Electronic multicriticality in bilayer graphene 
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We map out the possible ordered states in bilayer graphene at the neutrality point by extending the 
previous renormalization group treatment of many-body instabilities to finite temperature, trigonal 
warping and externally applied perpendicular electric field. We were able to analytically determine 
all outcomes of the RG flow equations for the nine four-fermion coupling constants. While the full 
phase diagram exhibits a rich structure, we confirm that when forward scattering dominates, the 
only ordering tendency with divergent susceptibility at finite temperature is the nematic. At finite 
temperature, this result is stable with respect to small back and layer imbalance scattering; further 
increasing their strength leads to the layer antiferromagnet. We also determine conditions for other 
ordered states to appear and compare our results to the special cases of attractive and repulsive 
Hubbard models where exact results are available. 



I. INTRODUCTION 

Understanding itinerant electronic systems with com- 
peting ordering tendencies is among the most profound 
challenges in today's condensed matter theory. In one- 
dimensional systems, powerful theoretical tools are avail- 
able for answering some of the questions^, but extending 
the techniques to higher dimensions has met with limited 
success. Often the problem is how to treat the various 
ordering tendencies on equal footing without an inherent 
bias towards any one of the possible ordered states. 

In this regard, bilayer graphene at, and near, the neu- 
trality point can be regarded as a model system. To a 
first approximation, there is a conduction band and a 
valence band that touch quadratically near two points, 
K and K' = — K, in the Brillouin zone 2 -! 3 -. Even when 
all electron-electron interactions are ignored, such a sys- 
tem would have low-temperature susceptibilities which 
diverge as ~ In T towards a number of different ordered 
states. While there are no known exact solutions, such 
a situation is expected to lead to instabilities with re- 
spect to infinitesimal electron-electron interactions. The 
challenge is then to identify the conditions under which 
any one combination of the various possible states gets 
preferably selected as the temperature is lowered. 

Since a many-body ordering appears already at weak 
coupling this problem is amenable to the renormaliza- 
ton group (RG) approach, whose advantage is that it 
can account for the competing tendencies in an unbiased 
way. Moreover, since the few-milli-clcctron-volt energy 
scales associated with ordering extracted from present- 
day experiments 4 ^— are much smaller than the natural 
upper cutoff in the problem originating from the split- 
off bands derived from the dimerized sites (~ 200 — 300 
meV), the physical system itself is expected to be well 
described by a weak coupling theory. Therefore, we ex- 
pect that the competition among the number of inher- 
ently strong-coupling phases can be accessed within such 
a weak coupling approximation. 

The previous RG treatments of this problem presented 
in Refs. HTIIT2I consisted of first building a low-energy ef- 



fective field theory, which, when electron-electron inter- 
actions are neglected, can be thought of as a Gaussian 
fixed point of the RG scale transformation^ 3 , with dy- 
namical critical exponent z = 2. Except under some 
non-generic fine-tuned initial conditions, contact interac- 
tions have been shown to be marginally relevant at this 
fixed point. Such four-fermion terms in the low-energy 
effective field theory arise from microscopic electron- 
electron interactions V ee (r) whose Fourier transform is 
non-divergent in the small wavevector limit. They could, 
for instance, correspond to 1/r Coulomb interactions 
screened by proximity to metallic gates. Within this ap- 
proach, the electronic modes with momenta in a thin shell 
(1 — A£)A < |k| < A near the cutoff A, and arbitrary fre- 
quency lu, are integrated out, while the change in the ef- 
fective action is monitored as the process is iterate d 13 ' 14 . 
To determine the leading instability, infinitesimal sym- 
metry breaking source terms were introduce d 11 ' 12 ' 15 and 
included in the process of renormalization. The source 
term with the strongest divergence was then identified 
as the most dominant ordering tendency. In the case of 
purely forward scattering, or, in the notation of this pa- 
per, for gA lg only, the leading instability was found to be 
toward the electronic nematic state. This state is gaplcss, 
with either two or four Dirac points near each K-point 
depending on the strength of the order parameter. In 
the case of the Hubbard model, there is additional back 
scattering, gE K = ^9A lg , and layer imbalance scatter- 
ing gA 2u = 9A lg 1 and the leading instability is found to 
be toward the layer antifcrromagnctic state. The single- 
particle (electronic) spectrum of this phase is gapped. 

In a similar approac h 16 ' 17 , the 1/r Coulomb inter- 
actions among the electrons in bilayer graphene were 
first screened using RPA, and then the full q- and in- 
dependent effective interaction was used as the initial 
condition for the subsequent Wilson-like RG treatment. 
While said approach can be criticized on the grounds that 
the screening of the long-range tail of the Coulomb in- 
teraction originates from integrating out electrons all the 
way down to the Fermi energy, which are double counted 
when reintroduced for the RG treatment, the results ob- 
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taincd using this approach are in qualitative agreement 
with the results obtained previousl y 11 ' 12 . 

To this end, we present an extension of the previous RG 
treatment of the problem to finite temperature 1 ^— and 
finite externally applied perpendicular electric field. This 
allows us to include the competition between broken- 
symmetry phases with gaps in the electronic spectrum, 
which may be energetically favorable, and gapless states, 
which may be entropically favorable. We also study the 
gradual suppression of an ordered state as the externally 
applied electric field is increased. Since temperature is 
treated explicitly, we can obtain the transition tempera- 
ture directly, without making any of the ad hoc assump- 
tions inherent in translating the value of the RG scale I 
at which the couplings diverge into temperature. 

As has been noted early on in the context of one- 
dimensional electron system o 21 ' 22 , it is very useful to 
compare the results of an approximate RG approach to 
known exact results 2 ^— . Despite the scarcity of exact 
results in higher dimensions, we can compare our re- 
sults to some of the non-trivial properties of the Hub- 
bard model at half-filling, which can be either estab- 
lished exactly 2 ^— or can be obtained from Monte Carlo 
simulations 3 - . - — . In this regard, it was shown in Ref. 
IT2I that starting with a repulsive Hubbard model on a 
honeycomb bilayer lattice at half-filling for U <C t± < t 
leads to the layer antiferromagnet as the most dominant 
instability. In this work, we confirm the previous finding 
using the finitc-tcmpcraturc RG scheme. We further es- 
tablish that, if we fix the value of the nine four-fermion 
coupling constants to correspond to the values derivable 
from the Hubbard model, then the low-energy effective 
field theory possesses the SO (4) symmetry of the micro- 
scopic Hubbard model 2 ^. As a consequence, for an at- 
tractive Hubbard model the result of our (approximate) 
RG analysis recovers the exact result that the s-wave su- 
perconducting order parameter can be continuously ro- 
tated to the "CDW" order paramete r 26 ! 27 . Since, for bi- 
layer graphene, the charge-ordered state does not break 
the discrete translational symmetry of the lattice, it is 
not strictly a density wave, but rather corresponds to 
the layer-polarized state (LP). This can be seen in our 
RG equations; the LP and s++ superconducting source 
terms are identical provided that we start with the val- 
ues of the four-fermion coupling constants corresponding 
to the Hubbard model with U -C t± < t. Moreover, 
since, in the weak coupling limit, we can map the micro- 
scopic lattice interactions to the four-fermion coupling 
constants in the continuum effective field theory, we can 
ask what happens when we add a &1-62 interaction V in 
addition to the on-site attraction U (see Fig. [1]) . When 
V is repulsive (attractive) we find that the exact degen- 
eracy between the LP and s ++ SC states is lifted in favor 
of the LP (s ++ SC) as expectec^ 3 .. 

Among the differences between our present approach 
and the related weak coupling approach employed in 
Refs. [l6lfl7l is the fact that we perform our analysis at 
finite temperature, which leads to different RG equations 



for the couplings than at zero temperature. In addition 
to the advantages mentioned above, this allows us to sys- 
tematically determine all possible outcomes of the RG 
equations in the nine-dimensional space of initial cou- 
plings. We also avoid screening the Coulomb interaction 
with the bilayer graphene low-energy degrees of freedom 
that enter our Wilson RG analysis. Rather, we assume 
that it is screened due to either finite temperature or the 
presence of external metallic gates. Finally, we do not 
rely on mean-field theory to determine the phases either 
directly from the bare couplings (i.e., without RG)2i~— , 
or on a rcnormalized mean field treatment^ 7 -. The short- 
comings of other approache a 34 ' 36 have been discussed in 
Ref. M- 

Within this formulation, as shown later in the text, the 
flow equations for the nin o 12 ' 16 coupling constants con- 
tain additional thermal factors, with an effective temper- 
ature T that grows under RG as e 21 . The flow equations 
()19l) for the coupling constants describe two competing 
tendencies — the term proportional to a product of two 
four-fermion coupling constants tends to enhance their 
growth, while the thermal factors suppress the flow of 
the coupling. For any fixed initial couplings and at a 
high enough temperature, the couplings saturate to fi- 
nite values as I — > 00. As the temperature is lowered, 
the coupling constants saturate at higher, but still finite, 
values. At the transition temperature, T c , the coupling 
constants diverge as £ — > 00. Below T c , the coupling 
constants diverge at a finite value of t. The effects of 
trigonal warping, parametrized by a velocity 1)3, can be 
readily included within this formalism as wel l 16 ' 17 . Like 
temperature, trigonal warping tends to suppress the flow 
of the couplings. As a result, even at T — 0, a critical 
coupling strength must be exceeded for a phase transi- 
tion to occu r 11 ' 16 . The strength of the critical coupling 
vanishes as V3 vanishes. 

The RG flow equations of the (infinitesimal) source 
terms for a multitude of symmetry-breaking order pa- 
rameters reveal that, at the transition temperature, the 
source terms A acquire an anomalous dimension, r/^. 
Analysis of the free energy correction to O (A 2 ) fur- 
ther reveals that, within this approximation, the phys- 
ical susceptibility for a particular A diverges as T — > T c 
if 77 a > 1. Using this condition, we determine the phase 
diagram for different initial couplings (see Fig. [5]). 

We find that, for purely forward electron-electron scat- 
tering, <7Ai„) the only order parameter with a divergent 
susceptibility at finite T is the nematic. Moreover, this 
is stable with respect to the presence of small, but finite, 
back scattering and layer imbalance scattering (i.e., 
the difference between intra- and interlaycr scattering) 
gA 2n ■ Performing the analysis at finite temperature is 
crucial for revealing this stability. Upon increasing the 
back and layer imbalance scattering, the only other di- 
vergent susceptibility is toward a layer antiferromagnetic 
(AF) state. Reversing the sign of the back scattering 
while fixing the layer imbalance scattering results in a 
quantum spin Hall state (QSH). Reversing the sign of 
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the layer imbalance scattering while fixing the sign of 
the back scattering gives us a layer-polarized state (LP). 
Reversing the sign of both may lead to an s-wave super- 
conductor. For small gA 2u /9E K and gE K ~ 9A lg > 0, we 
may find a Kekule current state (KC). These results are 
summarized in Fig. [5j which shows the phase diagram in 
the space of initial gA lg , 9e k , and gA 2u - 

Remarkably, the flow equations for the nine coupling 
constants can be analyzed in their entirety at T c . We 
find that, if a coupling constant diverges, it grows as 
e 2e . At the same time, the ratios of the coupling con- 
stants may either approach values determined by a two- 
parameter family of functions, which we call the target 
plane, or four isolated fixed ratios that do not belong to 
the fixed plane. For each of these cases we determine the 
symmetry-breaking channels with divergent susceptibili- 
ties at T c . The results are summarized in Fig. |H1 

The rest of the paper is organized as follows. In Sec. 
II, we present our model for the system. Section III is 
dedicated to the thorough analysis of the RG equations 
and our main results. Section IV deals with the effects 
of an applied perpendicular electric field on the phase 
boundaries. Our conclusions arc presented in Sec. V. We 
give details of our derivations in the appendixes. 



II. HAMILTONIAN 

We will be employing a low-energy effective theory 
for the bilayer graphene lattice. This lattice and the 
associated Brillouin zone are shown in Fig. [TJ Our 
model includes the nearest-neighbor intralayer hopping 

70 = t, the hopping between dimerized sites 71 = tj_, and 
the nearest-neighbor intcrlayer hopping between 11011- 
dimerized sites 73. It is this last hopping that is respon- 
sible for trigonal warping. Experimentally^, 70 ~ 3 eV, 

71 w 0.4 eV, and 73 sa 0.3 eV. Throughout this paper, 
we will use units in which Ub = ft = 1 ■ 




FIG. 1: (a) The honeycomb bilayer lattice formed by bilayer 
graphene. We represent the bottom layer, 1, with black circles 
and the top layer, 2, with red squares. The a; sites are the 
dimerized sites, and the hi sites are the non-dimerized sites. 
We include the nearest-neighbor intralayer hopping 70, the 
hopping between dimerized sites 71, and the nearest-neighbor 
interlayer hopping between non-dimerized sites 73. (b) The 
Brillouin zone associated with the honeycomb bilayer with 
the parabolic degeneracy points K = ^JL x and K' = K 
marked. 



a,\ site with a nearest-neighbor b\ site. The possible val- 
ues of S are — ^ax+ |ay, yax+ |ay, and —ay, where 
the lattice constant a rj 1.4 A. Whenever there is a sum 
on 5, we sum over these three values; if S appears without 
a summation over it, on the other hand, then we choose 
one of these three values for it. 

We may derive our low-energy effective theory for the 
above system by either projecting out the high-energy 
modes^ or equivalently by writing the above theory as a 
coherent-state path integral, integrating out the dimer- 
ized site o 12 ' 15 , and expanding around the K and K' 
points. The resulting theory is 

H = Ha + Hint (2) 



A. Non-interacting Hamiltonian 



where 



The tight-binding model for the lattice described above 



H, 



tb 



70 



J2 (at(R)MR + 8) + <4(R-)MR - 8) 



R,5.i 



+ h.C.) 

- liJ2(a{„(R)a 2<T (R) + h.c.) 

R,er 



73 (bL(n + 6)b 2a (R + 5 + S')+h.c), (1) 



R,5' 



where (a, 6) m(T (r) annihilates an electron on the (a, 6) 
sublattice on layer m and site r with spin a. The vectors 
R, are the positions of the dimerized sites within a unit 
cell, and S represents one of three vectors connecting an 



(3) 



|k|<A<r=T.4- 



In the above Hamiltonian, the Fermi spinor, which de- 
scribes the modes in the vicinity of the ±K points and 
concentrated at b sites in layers 1 and 2, is 



/ V&>(k) \ 



(4) 



The first of the two matrices in © describes the parabolic 
dispersion and the second is a linear term that results in 
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trigonal warping: 
# k = 

Hi 



find that 



1 



'k 

Jtw) 



2m . (« 

v 3 (k x A x 



kyky) 



QkxkylEjyj 



where 



A, 



lax, 2J S = T3CT2 



T3CT1, A y 



-lcr 2 - 



(5) 
(6) 
(7) 



(8) 
(9) 



In terms of the tight-binding parameters in our lattice 
Hamiltonian, the effective mass m* is 



and the trigonal warping velocity v 3 is 
v 3 = 3a7 3 . 



(10) 



(11) 



Experimentally^, m* « 0.029m e , while the value that 
we obtain from the above formula and the experimental 
values of the hopping parameters given above is m* m 
0.038m e . The value of the trigonal warping velocity used 
in fitting the experimental data£ is v 3 ps 1.41 x 10 5 m/s, 
while that obtained from the above formula is v 3 ps 1.91 x 
10 5 m/s (reference [HI assumes a value of v 3 — 10 5 m/s). 
The origin of these admittedly unimportant and small 
discrepancies is unclear at this time. 



B. Symmetry classification 

The space group symmetry operations which leave the 
Hamiltonian invariant at the T point form a point group 
Z?3d. Similarly, at the ±K points, the symmetry opera- 
tions form a point group D 3 . The character tables^ of 
these two groups are shown below. 



D 3d 


E 


2C 3 


3C 2 


i 


2S e 


3cr d 




1 


1 


1 


1 


1 


1 


M, 


1 


1 


-1 


1 


1 


-1 


E g 


2 


-1 





2 


-1 





Ai n 


1 


1 


1 


-1 


-1 


-1 


A 2u 


1 


1 


-1 


-1 


-1 


1 


E u 


2 


-1 





-2 


1 






D 3 


E 


2C 3 


3C 2 


Ax 


1 


1 


1 


A 2 


1 


1 


-1 


E 


2 


-1 






The sixteen 4x4 matrices that operate in the layer and 
±K valley space, can be grouped based on their trans- 
formation properties under these group operations. We 



Alg + 


1 4 






E g + 


(lo-i, r 3 o- 2 ) 


Axu- 


r 3 l 


A 2u + 


lcr 3 


E u — 


(T30-i,-lt7 2 ) 


AiK + 


T 1 a 1 ;T 2 (Jx 


A 2K - 


Tx0 2 ,T 2 a 2 


E K + 


(Til, -T 2 cr3; -T2l, -T1O3 



The ± next to the name of the representation denotes 
whether the particular operator is even or odd under time 
reversal symmetry. An equivalent classification can be 
found in Ref. [I?], though the notation is different. 



C. Interaction Hamiltonian 

As shown previously, assuming the microscopic lattice 
interactions are falling off faster than 1/r 2 , there are 
nine marginal interaction couplings at the Gaussian fixed 
point when T = and 1)3 = 0. The interaction term in 
the Hamiltonian is therefore 



Hi 



1 



E 

s 



9S_ 
2 



E E'^ka^k+q.aV'kV'^k' 



k,k' 



(12) 



The sum over S includes the 16 matrices belonging 
to the 9 representations. Since the couplings for the 
squares of the operators belonging to the same repre- 
sentation must be the same, we have 9 independent 
couplings. So, for example, for the E g representa- 
tion, the corresponding interaction term is schematically 

\dE g ((V^l^i^) 2 + (ipiT 3 a 2 i> a ) 2 y 

We may think of gE K as representing back scattering, 
9Ai g + <7a 2 „ as representing intralayer scattering, and 
9A lg — 9A 2u as representing interlayer scattering, as is 
demonstrated in our previous worki^. If we introduce 
a density-density interaction V(r) into the microscopic 
tight-binding Hamiltonian, then the forms of these cou- 
plings are given by Equations (119)-(121) in Ref. [TH, 
where they are denoted by <?Ai , 9d > and gp , respectively. 

As we will see shortly, if we start with these three cou- 
plings, then the other six will be generated under RG. 
All nine couplings may also be thought of as interactions 
between local fluctuations of different order parameters. 



III. FINITE-TEMPERATURE 
RENORMALIZATION GROUP 

We are interested in introducing the temperature T 
and the trigonal warping velocity v 3 explicitly into our 
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rcnormalization group transformations. To proceed, we 
rewrite the partition function as a coherent-state Grass- 
man path integral: 



To one-loop order, the RG flows of the coupling con- 
stants have the form 



(13) 



where 



So 



E 

n=-oo |k|<A o-=t,4- 



n is an integer, and the Matsubara frequency is lu u 
(2n + l)wT. The interaction term is 



(14) 



and 



1 1 

Mr,r) = -s E 7 E e"™" T e ik ' r ^K). (15) 



B ^ L 

n=-oo |k|<A 



Equivalently, we may write the interaction term as 
" ^° ^ 3=1 m=l 



(16) 



Note the absence of explicit spin subscripts on the (eight- 
component) ip = {ip^,if}\) T fields. The rj m ' ) matrices are 
defined in Eqs. (|C7)) - ([CT5|) . and rrij is the multiplicity of 
the jth representation. 

Our renormalization group procedure consists of split- 
ting the ip fields into fast and slow modes and progres- 
sively integrating out the fast modes with momenta re- 
stricted to the small shell A(l — A£) < |k| < A with no 
restriction on the Matsubara frequencies cj n . After each 
such mode elimination, we choose to rescale the momenta 
in the effective action for the slow modes such that the 

(2) 

new cutoff is again A and that the term is left in- 
variant. If we also wish to keep the iu n term invariant, 
and take A£ to be infinitesimal, we find that the temper- 
ature and the trigonal warping velocity flow under RG 

as 



dT 



dv 3 
dl 



2T =S> T{£) = e M T 
v-3 v 3 (£) = e e v 3 . 



(17) 
(18) 



In general, these flow equations will be corrected once 
interactions are taken into account, but for the couplings 
of choice here, the corrections appear only at two-loop 
order. 



dgi 



9 9 4 
3=1 k=l a=l 



where i, like j and k, extends over the aforementioned 
nine irreducible representations of the groups of the 
wavevector V and ±K; A^ k are listed in Appendix C. 
The dimensionless parameters that enter as the argu- 
ments of the $ functions are 



t(£) 



v 3 (i) 
A /2m* 

T{£) 
A 2 /2m 



(20) 
(21) 



The $ functions are determined by the integrals in Eqs. 
(154 HAT2)! . 

As shown in Appendix [XI these integrals can be evalu- 
ated explicitly when v 3 = in terms of elementary func- 
tions or when t = in terms of complete elliptic integrals. 

In the discussion that follows, we will make use of the 
asymptotic behavior in the limit of £ — > oo: 



-2( 



It 

T2t3 



+ ... for a = 1,2, (22) 
. for a = 3,4, (23) 



where t = t(0) is the initial dimensionless temperature 
and ". . ." represent terms that are smaller than the lead- 
ing terms. 



A. General analysis of the RG flows 



In general, the flow equations (fl9|) describe two com- 
peting tendencies. The term proportional to gjgk tends 
to cause an increase of the absolute value of the coupling 
constants as £ increases, while the $ functions tend to 
zero as £ increases due to the increase of their arguments 
y 3 {£) and t(£). Numerical analysis of the flow equations 
reveals that, for fixed values of the initial couplings and 
for a sufficiently large value of the initial temperature t, 
there is a certain value of £ where the flow becomes stag- 
nant and the coupling constants g tend to finite values as 
£ — >• oo. Therefore, if the initial couplings arc small, they 
remain small as long as the initial temperature is suffi- 
ciently large even as all the modes are integrated out . In 
this regime, weak-coupling RG is entirely justified. Low- 
ering the initial temperature, while keeping the initial 
couplings fixed, causes an increase of the value of the RG 
parameter £ where the coupling constants stop flowing 
and an increase in the limiting value of the coupling con- 
stants. At a critical initial temperature t c , the coupling 
constants g diverge as t oo. For an initial temperature 
t < t c , the coupling constants diverge at finite I. 



6 



The role of trigonal warping is to cause additional sup- 
pression of the increase of the absolute value of the cou- 
pling constants. Thus, for fixed initial values of the cou- 
pling constants and for sufficiently large initial 1)3, the g's 
do not diverge even at t = 0. 

Therefore, as stated previously^, for fixed initial 1)3, a 
critical value of the initial coupling(s) must be exceeded 
for a runaway flow of the coupling constant(s), which we 
associate with a phase transition, to occur. 

In order to understand the nature of the possible or- 
dering tendencies, we first analyze the asymptotic behav- 
ior of the equations (fT9")l when t = t c > and I — > 00. 
Provided that at least one coupling g r diverges, we have 
managed to enumerate all possible solutions for the sta- 
ble "rays" along which ratios with the other couplings 
9jldr tend to constants. The detailed analysis of these 
solutions is given in Sec. lIII Gl Along such a stable ray, all 
nine differential equations "collapse" onto one, namely, 



dg r t 2 e 



-2C 



H =A ^ 2t r 



■AS 



(24) 



Here, and in the remainder of the paper, if an index is in 
parentheses [e.g. (r)], then there is no automatic summa- 
tion over r unless explicitly stated. The coefficient *4( r ) 
depends on the stable ray along which the couplings di- 
verge and ". . ." denotes terms that vanish faster than 
e~ 2£ . Combining the asymptotic behavior of the $ func- 
tions as t -> 00, given by Eqs. (|22"]). ([2"3"]l. and ([15]). the 
coefficient may be expressed as 



9 9 



(25) 



j=l k=l a=l 



(r) 

where the = gj/g r is the ratio of two couplings along 
the stable ray. The solution of differential equation (|2"4")l 
is 



gr (i) = 



A 



, as £ —t 00 , 



(r) 



(26) 



where ". . ." denotes terms that are smaller than e 2i as 

£ ->• 00. 



B. Susceptibilities and the nature of the symmetry 
breaking 

To find out what symmetry-breaking tendencies dom- 
inate, we start by introducing source terms into our ac- 
tion: 

32 00 
AS = $>f A Y, E^K)0 (i VkK) + 



z—l n— — oo k 



16 



lE^ E E^K)o«c k (— n) 



+ C.C. 

(27) 



We may think of these terms as "forces" that couple 
to various observables, which acquire nonzero averages 
whenever the system enters the appropriate phase. Note 
that only 18 of the 32 particle-hole source terms intro- 
duced are symmetry inequivalent. Similarly, only nine 
of the 16 particle-particle source terms are inequivalent. 
The transformation properties of the former under the 
various symmetry group operations are summarized in 
Table |U Again, note the absence of explicit spin sub- 
scripts on the (eight-component) ip = (ip^,ipi) T fields. 
Terms such as i/j'Oip* should be understood as matrix 
multiplication, i.e., ^„ « =1 ^0^^. We will see later 
that only two of the particle-particle, or superconduct- 
ing, orders can appear, namely the A\ g and A2 U orders. 
These correspond to s++- and s_| -wave superconduct- 
ing orders, respectively. Both arc s-wave, but the s ++ 
order parameter has the same sign on both layers, while 
the order has opposite signs on each layer. To one- 
loop order, we find 

3=1 0=1 

j=l a=l 

where the (32 x 9) matrix and the (16 x 9) matrix 



din A? 
d£~ 

dlnA ? 
dT~ 



B\y are defined by Eqs. (|C29| - (^33l and (IC35l - (|C36| . 
Note that Eqs. (|28| and (|29[) can be readily integrated, 
and we find that 

Af /pp (£) = Af /pp (0)e 2£ cxp[ttf /pp (£)], (30) 

where 

«fW = EE4° fdl'g^aWW)), 

3=1 a=l J ° 

(31) 
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nrw = EE 



B 



(a) 



= 1 a=l 



d£' gj (£')$ a (v 3 (f),t(£')) 



(32) 

At t = t c > 0, as £ —> 00 the e 2e increase of a divergent 
coupling g r exactly balances the e~ 21 decrease of the <3> 
functions and the right hand sides of the above equations 
tend to constants, 



dlnAf 



dlnA pp 



= 2 



2B ph 



as 



00, 



as £ — > 00. 



(33) 



In other words, the engineering dimensions of the source 
terms, which are equal to 2, are corrected by the anoma- 
lous dimensions 

4 h/vp = —p— (34) 



A, 
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Group rep. 


Matrices 


Trans. 


TRS 


Inv. 


Mirror refl. (ad) 




Ai g charge 


1 4 eg 1 


+ 


+ 


e 


e 


Charge instability 


A2 g charge 


T303 ® 1 


+ 


- 


e 





Anomalous quantum Hall 37 ' 43 


Eg charge 


(1(71, T3CT2) Cg) 1 


+ 


+ 


e 


e/o 


Nematioii^ 


Ai u charge 


T 3 1 ® 1 


+ 


- 








Loop current 


A2u charge 


1ct 3 <g 1 


+ 


+ 





e 


Layer-polarized 34 1 36 


E u charge 


(T3<ri, -1<t 2 ) ® 1 


+ 


- 





o/e 


Loop current II (ME2) 


-A ik (Aig/Aiu) charge 


Ticri Cg) 1; T2CT1 Cg) 1 


— 


+ 


e/o 


e/o 


Kekule 4 ^ 


A 2 k {Aiu 1 Aig) charge 


Ti<72 Cg 1; T2CT2 Cg) 1 


— 


— 


o/e 


e/o 


Kekule current 


Ek (E g /E u ) charge 


(Til, -T2CT3) ® 1; (-T2I, -T1CT3) ® 1 


— 


+ 


e/o 


(e/o)/(o/e) 


Charge density wave 


Ai g spin 


I4 Cg (7 


+ 


- 


e 


e 


Ferromagnetic 


A 2g spin 


T3CT3 Cg) CT 


+ 


+ 


e 





Quantum spin HaUi^^ 


Eg spin 


(1cti,T3CT 2 ) <g>CT 


+ 


- 


e 


e/o 


Spin nematic 


Aiu spin 


T3I Cg ct 


+ 


+ 








Staggered spin current 


Aiu spin 


1(T3 Cg CT 


+ 


- 





e 


Layer AF^iML 


Eu spin 


(T3CT1, — 1(72 ) Cg a 


+ 


+ 





o/e 


Loop spin current II 


Aik (Aig/ 'Aiu) spin 


ticti Cg ct; T2(Ti Cg ct 






e/o 


e/o 


Spin Kekule 


A 2K (A 2u /A 2g ) spin 


Tia 2 Cg (?; T2CT2 Cg CT 




+ 


o/e 


e/o 


Spin Kekule current 


E K (Eg/Eu) spin 


(Til, T2 CT3 ) Cg CT; (— T2I, Tl CT3 ) Cg CT 






e/o 


(e/o)/(o/e) 


Spin density wave 



TABLE I: Table of all particle-hole phases considered, listed according to what representation of the D^d point group they 
transform. The Kekule and density waves have a wave vector of K. 



due to the electron-electron interactions. Again, in the 
above equation, there is no summation over r, which cor- 
responds to the divergent coupling g r that we divided by. 
The values of the i3's are 



k=l 



(35) 



(36) 



fc=i 



where B 1 - 1 /^ is given by the sum of Eqs. (jC30)) and (jC32|) 
and i^ 1 / 2 ' is given by (|C35|) . Note that the expressions 
for A( r ) and B p ^ pp depend on the choice of g r , but the 

In order to calculate the physical susceptibility toward 
various ordering tendencies, we calculate the correction 
to the free energy due to the presence of the symmetry 
breaking source terms^. We find that 



*/(A) 



(37) 



16 .00 4 

£ / d£e-^\Af(e)\ 2 Y,<^a(M£)Am 

„-_1 JO __1 



m 
16^ 



The a coefficients are given in Appendix IDlbv Eqs. (|D1|I - 
(ID41) . 

The susceptibilities are then simply given by second 
derivatives of the free energy with respect to the bare 



values of the appropriate source terms, 



ph 



X, 



xf = 



d 2 f 



d[A ph (£ = 0)} 2 '' 
d 2 f 



(38) 



d 2 f 



d[RcAf(^0)] 2 0pm Af (I = 0)] s 



(39) 



Using Eqs. (|30| and J37J), we find that the susceptibilities 
given above may be written as 



X, 



971 \ " ph 

= 8^< 



d^^^M,^))- 



(40) 



.ph/pp. 



0), being auxil- 



Note that the source terms, A£ 
iary fields, do not appear. 

Any divergence in the susceptibilities has to come from 
the regions of large £ in Eq. (|37|) where the asymptotic 
expressions derived earlier hold. Therefore, since, for t = 
t c > 0, the asymptotic behavior of the 3> functions is 
e~ 2e , the condition for the divergence of a susceptibility 
in a particle-hole or particle-particle channel i is 



v Ph/ PP > L 



(41) 



Next, we will relate the anomalous dimensions of 
the source terms rj pfl ^ pp to the susceptibility exponents 

ph/pp 
li 
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C. Susceptibility exponents 

For t > t c , sufficiently close to t c the asymptotic be- 
havior of the coupling constants is still approximately 
described by Eq. If wc integrate it from £ to £, 

both of which are asymptotically large (and temperature 
independent), but not infinite, then we find 



1 



1 



A 



g r (£,t) g r (£o,t) 



4/ 



-2i\ 



(42) 



At t c we have l/g r (£o,t c ) = Ar^e 2l °/At c and we can 
write the above equation as 



9r(e,t) 



1 



1 



1 1\ (A 

t ~ t 



(43) 



Since £q is finite, g r (£o, t) is analytic in t at t c and can be 
expanded as 



d 

g r (la,t) « g r (£ ,t c ) + (t- t c )—g r (io,t) 

at 



where ". . ." represents terms of order (t — t c ) and higher. 
Therefore 



where 



1 



d 



1 



dtg r (£ ,t) 



Clearly, the susceptibility for a particular order diverges 
if the condition (|4"Tj) is satisfied. Note that only if 

rjf p = 2 do the susceptibility exponents acquire their 
mean-field values. This is in general not the case here, 
as will be elaborated on in the next section. 

It is also important to stress that these exponents 
are obtained within the one-loop approximation of the 
fermionic theory and are therefore not expected to be 
accurate. They are also not expected to be equal to 
the one-loop exponents obtained within an e-expansion 
of the corresponding bosonic theory, with the Landau 
functional for the ordering field. The ultimate critical 
behavior is determined by the universality class of such 
a bosonic theory. As an example, the finite tempera- 
ture phase transition into the nematic state belongs to 
the two dimensional three-state Potts modelii universal- 
ity class for which the exponent 7 = 13/9 (see Ref. H§). 
However, within our one-loop fermionic RG treatment, 7 
does not exceed 2/3. Nevertheless, the exponents calcu- 
lated within the present approximation give us important 
information about the physical character of the dominant 
ordering tendency, without any a priori bias toward any 
given order. 

Next, we will explicitly calculate the RG flows us- 
ing numerical integration of the RG equations for the 
couplings, the symmetry-breaking vertex terms and the 
physical susceptibilities. We do so for any interac- 
tion that can be written as J2qJ2i j Vij( c O ri i( c l) n j(~ c l)> 
as l _ » oo 1 t+(45) where Vy(q) is finite for any q; % and j run over sub- 
lattice and layer indices. To leading order in small V, 
q is either near or near ±2K. Such a microscopic lat- 
tice interaction will initially lead to only three of the 
nine four-fermion coupling constants in the low-energy 



A 



O) e -2l 



4/ 



, (44) 



•"V) -21q 

Atl 



Note that c r A( r ) > since g r (£o,t) increases in magni- 
tude as t ^ t+. 

The flow of the source terms at large £ at t > t c is de- 
termined by substituting the above result into the Eqs. 
(f2"5| - (f2"9"|) and taking the asymptotic limit of the $ func- 
tions at large I: 



din A?* 

d£ 

d£ 



= 2 



= 2 



B 



ph 
i(r) 



-21 



It 



-21 



It 



Cr(t-t c ) + ^e-™ 



(47) 
(48) 



Integrating from £q to £ and substituting to Eq. (|37| . we 
find that the singular contribution to the susceptibility 
for the symmetry breaking source term Aj is 



X 



ph/pp 



(t-t c y 



ph/ pp 



where 



ph/pp 



ph/pp 



(46) effective field theory being fmiteia, i.e., gA lg \ e=0 7^ 0, 

9A 2u \ e=0 ^ 0, gE K \ e=0 ^ 0. 

D. Forward scattering limit: nematic 

For Vij (±2K) = and equal inter- and intralayer in- 
teractions, the only non-zero bare interaction is gA lg ■ 
Without scattering between the K and K' valleys at the 
microscopic level, no new scattering between them can 
be generated in the RG flow, Eq. (flT)|) . In other words, 
9A 1K , 9a 2 ki an d gE K remain zero. The only couplings 
that appear in this case are those in the g and u rep- 
resentations. We studied the problem numerically and 
present the main results in Figs. [5] and [3] 

For any fixed initial U3, we find that there is a critical 
value of gA x below which the weak-coupling RG con- 
verges and no phase transition occurs, even at T = 0, as 
shown in Fig. [2^.. In this phase, where no symmetry is 
broken, there are four Dirac points in the vicinity of each 
K point, three of which are anisotropic and one, centered 
at K or K', which is isotropic. 

Above the critical value of 9a 1s , the only susceptibil- 
ity that diverges as t — >• i+ > is toward an electronic 
(50) nematic, i.e., toward a spin-singlct order parameter that 
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FIG. 2: (Color online) Phase boundaries for bilayer graphene 
with forward scattering only, (a) At finite trigonal warping, 
v-i, — 2m*vs/A and T = 0, the bare gA lg must exceed a 
critical value, given by the red line, in order for the system 
to enter a broken symmetry phase. Along the dashed line, 
vs, — 0.178, which is the value used to fit the experimental 
data in Ref. @. (b) The transition temperature as a function 
of the initial value of gA lg at 1/3 — 0.178. At any finite t c , the 
only susceptibility that diverges corresponds to the nematic 
order parameter (E g charge) , as shown in Fig. f4] 



transforms according to the E g representation (sec Table 
H]). We therefore conclude that, immediately below this 
temperature, the system enters this symmetry-breaking 
phase. In Ref. the experimental data is fitted using a 
value of the trigonal warping velocity corresponding to 
our dimensionless parameter 1/3 = 0.178. Here, and in 
the remainder of the paper, we use this value. The phase 
boundary for that particular choice of 1/3 is shown in Fig. 
[2jb). In Fig. [4] wc show the susceptibilities as a function 
of temperature for various order parameters. This plot 
corresponds to the vertical arrow in Fig. &h) where only 
the nematic susceptibility diverges. We find that the sus- 
ceptibility toward the nematic order parameter, despite 
being smaller than the others at large T, outgrows all 
competing susceptibilities as the temperature is lowered 
toward T c , and is the only susceptibility to diverge at T c . 



0.10 




4n /m 



OffiO 



FIG. 3: (Color online) The phase diagram for different val- 
ues of the forward scattering coupling, gA lg , trigonal warping 
velocity, vz — 2m*vz/A, and temperature t — 2m*T/A 2 . At 
any finite t c and at any value of 1/3 the nematic susceptibil- 
ity is the only one that diverges. The lighter cyan line, also 
shown in Fig. Ob), corresponds to vz = 0.178, which is the 
value used throughout the paper. 



While our analysis of the susceptibilities reveals that 
at t c > only the nematic susceptibility diverges for any 
fixed u 3 , this susceptibility does not diverge when ap- 
proaching the critical gA lg from below exactly at t = 0. 
Instead, we find that the susceptibilities for the layer an- 
tifcrromagnct (AF) and quantum spin Hall (QSH) order 
parameters diverge with equal exponents. This suggests 
that, at < t < t c , the system orders into a nematic 
state, while at t = there may be a coexistence of this 
state with AF and/or QSH^. 

The complete phase diagram for different values of 
ffAi j ^3, and t c is shown in Fig. [3j For the entire range 
of ^3's shown, we always find the nematic as the leading 
instability at finite temperature. 

In order to facilitate the comparison with the previous 
work, which also deals with the forward scattering limit 
at T = 773 = 0, we first note that the three couplings g\, 
g 2 , and g 3 in Ref. El correspond to g Alg , g A2g , and g Eg , 
respectively. Because V3 = 0, none of the other nine four 
fermion couplings are generated under RG. The "suscep- 
tibilities" calculated therein correspond to the logarith- 
mic prefactors on the right-hand sides of Eqs. (15) and 
(16) in Ref. [n]> an d are analogous to the more general 
expressions in Eqs. (J2SJ) d2HJ) of this publication. The 
physical susceptibilites, considered in this paper, can be 
straightforwardly obtained from such expressions by sub- 
stituting the flow of the source terms into the formula 
([57]). At T = 0, the divergences appear at finite I = t* , 
which can be set as the upper limit on the integrals in 
Eq. (|57|) . The coupling constant ratioii gA lg /gE g —> 
as £ — ¥ £* . The ratio gA 2 I 9e can approach cither 
mi ~ —0.525 or 771 3 w 13.98, i.e., the minimal or the 
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maximal root of the cubic equation (x — 14)a; 2 + 4 = 0. 

The analysis of the physical susceptibilities for the con- 
ditions stated in Ref. fill reveals that, as I — s- £*, the only 
physical susceptibility that diverges when gA 2 j 3e — > 
mi is towards the nematic; the others remain finite at 
i* . Similarly, the only susceptibility that diverges when 
9a 2 /9e ;i W3 is toward the quantum anomalous Hall 
state (QAH). Very recently, Fan Zhang et al. posted a 
preprint^ in which they recovered the T = i>3 = flow 
equations for the three couplings in Ref. [ill Under equiv- 
alent conditions to those in Ref. [n], they claim to have 
calculated susceptibilities and that the "strongest diver- 
gence occurs for the flavor spin channel broken inversion 
symmetry." These results are at odds with our findings. 
We believe the discrepancy to originate from their Eq. 
(24), which does not properly account for renormaliza- 
tion of composite operators (see Refs. HUHF 
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FIG. 4: (Color online) Various susceptibilities calculated from 
the free energy given by Eq. (|37[) with forward scattering only. 
Although the nematic susceptibility is lower than the others at 
higher temperatures, it is the only susceptibility that diverges 
as the temperature is lowered towards t c > 0. Here, — 
0.178, the bare couplings are gA lg = 0.161 x 4ir/m*, and all 
others zero. In this case, t c = 0.01. 



E. General density-density interaction 

In the previous section, we have shown how a sys- 
tem with forward scattering only at the bare level or- 
ders into the nematic state at any finite temperature. In 
general, however, other four-fcrmion coupling constants 
may be non-zero as well. In a previous work on this 
subject^, two of us showed how to find the bare interac- 
tion strengths corresponding to a screened interaction in 
the weak coupling limit. In addition to gA lg , two other 
couplings, gA 2u and <7e k , appear at £ = 0. Due to the 
presence of these couplings, all (3 functions are non-zero 
and all nine couplings allowed by symmetry are generated 
in the RG flow. 

Instead of seeking a critical temperature for a given 
set of initial couplings, we invert the procedure by fixing 



the transition temperature to t c = 0.01. This value is in 
accordance with the experimentally observed symmetry- 
breaking energy scale of ~ 2 meV. We then determine 
what values of the initial couplings would make the RG 
flow divergent at this temperature. This set of points de- 
fines a two-dimensional surface in the three-dimensional 
space of initial {gA t , 9a 2u , 9e k )- For each point on this 
surface, we find the list of phases for which the suscep- 
tibility divergence criterion, Eq. (|4Tj) , is satisfied. For 
certain initial conditions it happens that two or more sus- 
ceptibilities diverge at t c . In such situations we list all 
the phases with diverging susceptibilities (e.g., "N+AF" 
represents the region of initial couplings where both xn 
and xaf diverge, although not necessarily with the same 
exponent). Because our formalism is valid only for t > t c , 
the resulting state may be either one of the listed phases 
or a coexistence of several of these phases. In order to 
decide which phase (s) is present, one needs to construct 
a theory valid below t Cl such as a Landau theory with 
multiple order parameters. This is beyond the scope of 
the present paper. 

The phase diagram we find is presented in Fig. [5] 
One should understand the axes on this plot as fol- 
lows. When the microscopic interaction has a long range, 
the bare values of <?a 2 „ and gE K are negligible rela- 
tive to gAi ;l ■ They become larger for interactions with 
shorter ranged. For monotonically-decreasing repulsive 
interaction potentials, these two couplings do not ex- 
ceed gA lg and g^/2, respectively. gA 2u and gE K reach 
these values in the Hubbard limit, where the only micro- 
scopic interaction term is on-site. We therefore restrict 
our phase diagram to positive initial values of gA lg , to 
\9E K \/gA lg < 1/2, and to \gA 2u \/9A lg < 1- 

In the given range of initial couplings, we find a rich 
phase diagram with the following phases: 

(a) Nematic (N): This phase is stable for predomi- 
nantly forward scattering, i.e., when both gA 2u and gE K 
are small at the bare level. If one of these couplings re- 
mains small and the other becomes comparable to 9A lg 
the nematic susceptibility is still divergent, although 
other susceptibilities will diverge at these initial values 
as well. The nematic state is gapless, but it reconstructs 
the low-energy spectrum such that two out of four Dirac 
cones in each valley become gapped. 

(b) Layer antiferromagnet (AF): This phase occurs 
when all three bare couplings are repulsive and compara- 
ble, corresponding to a short-range repulsive interaction. 
In this state the magnetization on each undimerized site 
is finite, with the magnetization within one layer point- 
ing in one direction, and that in the other layer in the 
opposite direction. 

(c) Layer-polarized state (LP): This phase is preferred 
when the intcrlaycr repulsion is stronger than the in- 
tralayer repulsion (gA 2u (£ = 0) < 0), and the backscat- 
tcring is cither repulsive or weakly attractive {gE K (£ = 
0) > gA 2u {£ = 0)). In this phase, which is gapped, there 
is an imbalance of the electron occupation number be- 
tween the two layers. One layer is more occupied and 
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two plaquettes carry a circulating current, both in the 
same direction. This phase is gapped. 

For a graphical illustration of some of these phases, see 
Fig. 2 in Ref. [H 

In the entire plot the values of the bare dimensionless 
couplings m*gi/A-K for which the system orders at the 
preset t c = 0.01 and V3 = 0.178 are always smaller than 
0.15, which justifies our weak-coupling approach. 

The situation does not change qualitatively with varia- 
tions in temperature or in the absence of trigonal warping 
— we explored a range of temperatures 0.005 < t c < 0.02 
with and without trigonal warping and found that the 
general structure of the phase diagram in Fig. [5] does not 
change appreciably. This is not a coincidence and, later 
in the paper, we will map the phases that we may obtain 
by analyzing the behavior of the flow equations in the 
large £ limit, where trigonal warping is irrelevant. 



F. The Hubbard model — "hidden" symmetry 



FIG. 5: (Color online) The phase diagram of bilayer graphene 
with trigonal warping. The transition temperature is fixed 
to t c = 0.01 and 1/3 — 0.178. Predominantly forward scat- 
tering favors the nematic (N) phase. When backscattering 
and/or the difference between inter- and intralayer scattering 
is considerable at the bare level we find other phases: the 
layer antiferromagnet (AF), the layer-polarized state (LP), 
the quantum spin Hall state (QSH), the s++ superconduct- 
ing state (s++ SC), and the Kekule current state (KC). In re- 
gions where two or more susceptibilities diverge at the same 
t c we use "+" to denote a "coexistence" of multiple possi- 
ble phases. Whether the listed phases truly coexist or one of 
them is preferred must be determined from the full Landau 
function. Such an analysis is beyond the scope of this paper. 



the opposite layer is equally less occupied with respect 
to the symmetric, high-temperature, state. 

(d) Quantum spin Hall state (QSH): This state is pre- 
ferred when the backscattering is attractive (gE K {(■ = 
0) J$ 0), but if <7a 2 „ is attractive as well, it must be 



weaker (gE* 



0) < 9A 2 



0)). In this state, which 



is gapped, there is a spin current around each plaquettc 
circulating in the same direction in both layers. 

(e) s ++ superconductor SC): The conditions for 
this phase are similar to the previous one in that the 
backscattering must be attractive, but it must also be 
roughly stronger than (attractive or repulsive) |<M 2u I a t 
t = 0. This state opens a superconducting gap in both 
layers with the same sign on each layer. 

(f) Kekule current phase (KC): This phase appears in 
a thin sliver of initial couplings for which backscattering 
is repulsive and comparable to gA lg , while the inter- and 
intralayer couplings are roughly the same (<?a 2 „ = 0) ~ 
0). It breaks lattice translational symmetry and time- 
reversal symmetry. In this phase a supercell, consisting of 
three regular unit cells, is formed. Within the supercell, 



As an important check, we apply our RG procedure to 
a special case, namely the Hubbard model, about which 
we already know certain exact properties. At half filling, 
this model has a dynamical SO (4) symmetry^ on a bi- 
partite lattice. This symmetry is present regardless of the 
sign of U and the dimensionality When U is negative, 
this symmetry is particularly useful because the electrons 
have a tendency toward pairing. One good variational 
ground state for the negative U Hubbard model on a 
square lattice is a charge density wave, where one sub- 
lattice has a higher occupation number than the other. 
Another ground state that exhibits pairing is the s-wave 
superconductor. The pseudospin symmetry rotates be- 
tween these states. At half filling the tendency towards 
the charge density wave order must therefore be the same 
as the tendency towards the s-wave superconducting or- 
der. 

In the case of bilayer graphene, the nomenclature is 
slightly different. A difference in the number of elec- 
trons in one sublattice compared to the other corresponds 
to a layer-polarized state and not to a CDW because 
the layer-polarized state does not break the translational 
symmetry of the lattice. Among a plethora of super- 
conducting orders in bilayer graphene, the one that is 
obtained by pseudospin rotation from the layer-polarized 
state is the s ++ superconductor. Following the argument 
in the previous paragraph, we conclude that the tenden- 
cies towards the layer-polarized and s ++ superconducting 
orders are exactly the same in bilayer graphene at half 
filling with an attractive Hubbard interaction. 

In addition, a repulsive Hubbard model is related to 
its attractive counter part with an equally strong inter- 
action. The mapping between the two is given by22, 



c T (R) 



c T (R), 



(-l)M(R) 



(51) 
(52) 
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FIG. 6: (Color online) Susceptibilities for the repulsive Hub- 
bard model on the honeycomb bilayer. The AF susceptibility 
is the most dominant at all temperatures. The dominance 
of the AF susceptibility over the others is a sign of growing 
AF correlations. Here, vz = 0.178, the bare couplings are 
9A lg = 9a 2 „ = 3-E K = 0.0560 x 4n/m*, and all others zero. 
In this case, t c = 0.01. 



It is easy to demonstrate that, under these transforma- 
tions, the kinetic term of a Hubbard model on a bipartite 
lattice at half filling remains the same, while the interac- 
tion term changes sign. 

The pseudospin symmetry, as well as the mapping of 
the repulsive Hubbard model to its attractive counter- 
part, are also present in the honeycomb bilayer lattice. 
The question is then whether any of these properties sur- 
vive in the low-energy effective field theory, in which the 
only degrees of freedom considered are those around K 
and K'. In the j3 functions, Eq. (fH))) , the pseudospin 
symmetry is not apparent. Moreover, once we start from 
a set of bare interactions corresponding to the Hubbard 
interaction, all nine couplings are generated. However, as 
we will now demonstrate, both the pseudospin symmetry 
and the connection between the attractive and repulsive 
Hubbard models are present in the RG flow. These man- 
ifest themselves through certain linear combinations of 
the couplings that are invariant when the bare interac- 
tions correspond to the Hubbard model. 

We start by rewriting the mapping given by Eqs. (|51[) 
and ([5^1) for our low-energy effective theory. Fields with 
spin up transform as 



while those with spin down transform according to 
<,( w ) f„\ , J.( bl ) 



r(61/62). 



^(q) -+ <K J K-q), 



«i(q) 



«I(-q)- 



y TKi 



(54) 
(55) 



This mapping leaves the kinetic term, Eq. ([5]), invariant, 
as it should, but it changes the interaction term, Eq. (TT2"1) . 



\ ]T gs {^S^f -4^> (ftS$) 2 



(56) 



where each coupling constant gs is a linear combination 
of the coupling constants gs before the transformation. 
We find that four of the nine coupling constants, gA 2g , 
9E g , 9A lu , and g AlIll do not change sign, i.e., g L -> g i} 
under this mapping. We call these coupling constants 
"even." In addition, there are two linear combinations, 
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which change sign under the mapping, i.e., g a ^ — > —g a /b- 
We will refer to these as "odd." Finally, there are three 
remaining linearly independent combinations of the cou- 
pling constants, 



s 9i = 9A lg - 9A 2u ~ 2gE u , 

5 92 = gA lg - gA 2u - 25a 2K , 

593 = 9A lg +gA 2u -^gE K , 



(59) 
(60) 
(61) 



which arc neither "even" nor "odd," as they generate 
terms proportional to themselves as well as terms pro- 
portional to "even" and "odd" coupling constants in the 
RG flow. Clearly, a non-zero value of any 5gi would spoil 
the connection between the two models. 

In the Hubbard limit, we notice that all three Sgi's 
are zero at i = 0. We also see that, initially, all other 
couplings are exactly zero except for g a ~ 9U/4 7^ 0. 
Therefore, Eqs. (|5U)) - ([55]) . when applied to a repulsive 
Hubbard model with interaction strength U, changes the 
sign of the only non-zero coupling g a — — g a and leaves 
all the other couplings zero. This is precisely the bare 
interaction of an attractive Hubbard model with the same 
interaction strength. 

So far, we have shown that the low-energy effective 
field theories for the repulsive and attractive Hubbard 
models in bilayer graphene map onto each other, but only 
at the bare level. To show the equivalence at any £, we 
look at the flow of the coupling constants, i.e., their linear 
combinations, Eqs. (|5"71) - (|6Tj) . For the three couplings, 
Eqs. ([g9 ]) -([6l l . we find that 



-^r = Ptgi {5gi,8g2-,5g d ) 5g ^ 



0. 



(62) 



The last arrow means that all three ft functions vanish 
when all three Sg^s are simultaneously zero. Since this 
is true at I = in the Hubbard limit, it follows that no 
SgiS can be generated in the RG flow. 

On the other hand, the "even" and "odd" coupling 
constants do flow under RG, but there is a special struc- 
ture to their ft functions. The four "even" couplings flow 
according to 
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The flows of the two "odd" couplings are given by 



(o) 



d( 



(64) 



=1.7=1 fc=l 



The fact that all Sg^s are zero for any t has already 
been incorporated, and thus there are only six indepen- 
dent couplings. The structure of Eqs. and (fM)) 
makes them manifestly invariant under the transforma- 
tion, ([55)1 Therefore, an RG flow obtained within 
the continuum low-energy effective field theory corre- 
sponding to an attractive Hubbard model is described 
by the same set of differential equations as its repulsive 
counterpart. The flows of the "even" couplings arc identi- 
cal for the two cases, while those for the "odd" couplings 
differ only by a sign. The three couplings Sgi are all zero 
at any £ for both cases. Had the flow started from a point 
where at least one of the 5gi 's were finite, this correspon- 
dence would have been spoiled because additional terms 
appear in Eqs. (|63[) and (l64l) . 
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FIG. 7: Susceptibilities toward the layer-polarized and s++ 
superconducting orders in the attractive Hubbard model. In 
this case, both susceptibilities are equal. The bare coupling 
constants used here are gA lg = gA 2u — %9e k — —0.0560 x 
47r/m*, with all others zero. The transition temperature in 
this case is t c = 0.01. 

An immediate consequence of the mapping described 
here is that physical quantities obtained in our one-loop 
RG method for a repulsive Hubbard model are related to 
those obtained from its attractive counterpart. For exam- 
ple, the "critical temperatures" t c for the two models are 
the same. Of course, the layer antiferromagnetic phase 
for the repulsive Hubbard model and the layer-polarized 
and s++ superconducting phases for the attractive Hub- 
bard model both have the zero transition temperature 
because a continuous 0(3) symmetry cannot be broken 
at any finite temperature in two dimensions. The finite 
t c that we obtain within this approximate technique cor- 
responds to a gap scale that must be the same for both 
models. 

Having demonstrated these special properties of the 
RG flow in the Hubbard limit, we now compare the sus- 
ceptibilities for the layer-polarized and s++ supercon- 



ducting states. The a coefficients for the correspond- 
ing source terms in the free energy, Eq. (|37[) . are equal. 
Therefore, it is sufficient to look at the difference of the 
right hand sides of Eqs. and 
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Since none of the Sg^s are generated in the RG flow in the 
Hubbard limit, the source terms for the layer-polarized 
and s++ superconducting states flow in exactly the same 
way, their susceptibilities must diverge with the same 
exponent. This proves that our one- loop RG treatment 
respects the pseudospin symmetry of the Hubbard model 
at half filling. This argument remains valid for any value 
of the trigonal warping, which does not break particle- 
hole symmetry, since no assumptions were made about 
$ functions. 

Notice that the Hubbard limit is not the only case in 
which the mapping and consequent pseudospin symmetry 
are realized. Any model in which the bare values of all 
three Sg^s are simultaneously zero will exhibit the above 
correspondence as well. However, if we restrict ourselves 
to microscopic density-density interaction Hamiltonians, 
in which case only three of the four-fermion coupling con- 
stants arc initially non-zero, the pseudospin symmetry is 
present only if g Alg = g A2u = 2g Ell - 

At the end of this section we present numerically ob- 
tained susceptibilities for various orders in the case when 
gAi g = gA 2u = 2gE K - In Fig- El susceptibilities to vari- 
ous orders in the repulsive Hubbard model are shown as 
functions of temperature. The AF susceptibility domi- 
nates and is the only one to diverge at T c . Reiterating 
what was stated before, this divergence is an artifact of 
our one-loop RG approximation. Nevertheless, one-loop 
RG correctly singles out the state that is known to be 
the ground state at T = 0. The t c that we find should 
be thought of as a gap scale for the AF order. Figs. 
[7J and [5J compare the layer-polarized and s ++ supercon- 
ducting susceptibilities. In Fig. [JJ an attractive Hubbard 
model is studied. In Fig. [51 we consider the same model 
with an additional small 61-62 repulsion (left panel) or 
attraction (right panel). This additional term violates 
the pseudospin symmetry of the Hubbard Hamiltonian. 
In this case, the coupling constants are the same as in 
the Hubbard model, except that now <M 2u = j^pf gA lg , 
where e = V/U and V is the microscopic 61-62 interaction 
strength. When this interaction is repulsive, the system 
favors the layer-polarized state over the superconduc ting 
state analogous to, for example, the findings of Ref. l33l . 
Conversely, when this interaction is attractive, it favors 
derealization of the electron pairs and the concomitant 
superconducting ground state. Our numerical results are 
in the agreement this. 
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FIG. 8: (Color online) Comparison of the susceptibilities for the layer-polarized and s++ superconducting orders in the attractive 
Hubbard model with an additional small 61-62 interaction. In this case, <7A 2u = y^pf 9A lg , where e = V/U and V is the microscopic 
61-62 interaction strength. In both cases, the bare gA lg = ^9e k = —0.05 60 x 47r/m*. (Left panel) When a small 61-&2 repulsion 
is added, the layer-polarized state is preferred. Its susceptibility diverges at t c , while that of the superconducting state, as 
well as all other order parameters considered, reaches a finite value. Here, we take e = 0.1, in which case t c = 5.18 x 10 -3 . 
(Right panel) In the case of a small 61-62 attraction, the opposite is true — the susceptibility for the s++ superconducting state 
diverges, while that for the layer-polarized state remains finite at t c . Here, e = —0.1, in which case t c = 0.0288. 



G. Fixed ratios and broken symmetry phases 



become 



The list of phases found numerically in the previous 
section shows ordering trends for bilayer graphcne only 
for a certain kind of microscopic interaction that we be- 
lieve is relevant in a realistic system. The question is 
whether there are other possible ordered states in bilayer 
graphcne when all 9 symmetry-allowed couplings are in- 
cluded at I = 0. One way to answer this question is 
to numerically explore the entire 9-dimensional space of 
bare couplings for various trigonal warping parameters. 
Such an approach, although straightforward, would re- 
quire immense computational power and might even miss 
certain phases that are realized only for specific bare in- 
teractions. Fortunately, there is another approach to the 
problem that we discuss in this section. Instead of con- 
centrating on the bare interactions, we look at what hap- 
pens to the couplings and susceptibilities at large I. This 
allows us to enumerate all the possible phases regardless 
of the initial interactions. 

Previously, we discussed the asymptotic behavior of 
the RG equations at t = t c > 0. We know that at least 
one coupling will diverge as g r (() ~ e 2e - We divide all 
the other couplings by that particular coupling and find 
the (3 functions for the ratios, = gj/g r , to be 
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Here, a dot over a coupling constant represents a deriva- 
tive with respect to I. In the large £ limit, these equations 
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We now ask if these equations have any fixed points, or, 
in our terminology, "fixed rays" . These are obtained by 
demanding that the right hand sides of all 8 equations 
(|67p are simultaneously equal to zero. After finding the 
fixed rays, we need to determine whether each ray is sta- 
ble, unstable, or mixed by analyzing eigenvalues of the 
stability matrix Sjk = dpj /dpjp. Since Am is already 
defined in Eq. ([25)1 in terms of the ratios, the entire sta- 
bility matrix has well-defined eigenvalues for each "fixed 
ray" solution. In addition, the sign of Ai r ) determines 
the sign of the diverging coupling that we divide the oth- 
ers by; see Eq. 

If we find that a ray is stable, then, if we start with 
the coupling constants sufficiently close to the fixed ray, 
then the ratios of the couplings approach the given set 
of values as I — s- 00. Such a flow leads to a divergent 
susceptibility in at least one channel. If a ray is mixed or 
unstable, then, in the absence of fine-tuning, the RG flow 
cannot take the couplings toward such a ratio; even if the 
flow starts in such a direction for small £, it will be redi- 
rected toward some other ray that is stable. We therefore 
conclude that all the solutions that have even one pos- 
itive eigenvalue in their stability matrix are physically 
irrelevant. It is possible that some rays are marginal in 
certain directions, meaning that some of the eigenvalues 
of the stability matrix are zero, and stable in others. We 
do, in fact, find such physically relevant solutions. 

Following the procedure described above for all pos- 
sible choices of the divergent coupling, we find that the 
stable solutions of the RG flow are situated cither on a 
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FIG. 9: A plot of all of the phases found in the fixed plane described by Equations (170[) - (|72[) . We find nematic (N, E g charge), 
Kekule (K, A\k charge), spontaneous current, or magnetoelectric (ME2, E u charge), layer-polarized (LP, A2 U charge), Kekule 
current (KC, A2K charge), staggered spin current (SSC, A\ u spin), antiferromagnetic (AF, A2 U spin), quantum spin Hall (QSH, 

A2 g spin), s++ superconductor (s++ SC, A\ g singlet), and superconductor (s+_ SC, A2 U singlet) states. In addition to 

this fixed plane, we also find four isolated fixed points, which are described in the text. 



manifold that we call the "target plane" or on one of ratios: 
four isolated fixed rays. The "target plane" represents a 
set of stable rays that are marginal in two directions and 
stable in six others. The target plane and the phases cor- 
responding to each point within are shown in Fig. [51 We 
parametrize the plane in the following way. We choose 
as our parameters the following two coupling constant 
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Since, for certain fixed rays, gs u and/or gE K diverge, 
while gE g does not, these parameters take values in the 
interval (—00,00), including infinite values. With the 
chosen parametrization, we express each coupling at large 



16 



£ as 



9A lu 

9A 1K 



9M 



G{£) 



Ay 2 



C(x,y) 

Axy 
C{x~y) 



0, 9a 2 
G{£), 9A 2K 



(1 + x + 2y)< 
C(x,y) 

Ax 

t=t c (x,y) 



G(£), 9E q 



G(i), g Eu 



4y 



C(x,y) 



G(£), 



-2(l + x + 2y) 



9E* 



C{x,y) 




-2x{l + x + 


2y) ( 


C(x,y) 




-2y(l + x ■+ 





t=t c 



C(x,y) 



G{£), 
G(£), 
G(£), 



(70) 
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r 



where C(x,y) is a square root of a quartic polynomial 



E 



i9 3 



is a pos- 
oo. The 



and the "overall" coupling G(£) 

itive definite function of £ that diverges as £ 
expression for C(x,y) can be readily obtained from the 
definition of G{£), but is unwieldly, and thus we do not 
include it here. The ratios of any two couplings at large £ 
depend only on x and y, although sometimes these ratios 
may be infinite. 

In the special situation in which the parameters x and 
y are infinite, but their ratio is finite, we may reparam- 
eterize x and y as x = Rcosn and y = R sin rj and take 
the limit as R — > oo. The only diverging couplings in this 
case are gA 2g , 9A lu , 9e u , gA 1K1 9e k - Note that, for each 
rj, we obtain the same stable ray at rj + tt. Due to the 
fact that any two opposite points at infinity on the tar- 
get plane are identical, we conclude that the target plane 
is homeomorphic to a projective plane M.P 2 . In Fig. [51 
some of the phases, such as QSH, have hyperbolic phase 
boundaries and appear to exist in two disconnected parts 
of the phase diagram. However, due to the fact that the 
opposite points in the target plane are identical, these 
may be regarded as single and simply connected regions. 

The values of pj = gj/g r are readily obtained from 
([T0")) - ([T2")) . Without loss of generality we now set g r = gE c 
in Eq. ([26]). We obtain 



c 3 + 2x + 3x 2 + Ay + Axy + 8y 2 m* 

Aie ) = —6 . (73) 

( 9 > l + x + 2y Ait v 1 

We may obtain Bi,(E g ) from Eqs. ([33]) and ([3"6"j) . We 
can now determine the anomalous dimensions of the 
symmetry-breaking source terms defined in Eq. (|34[) . Re- 
markably, we see that the anomalous dimensions are con- 
tinuous functions of the two parameters x and y. For 
each point in the target plane, we determine the phases 
for which r\i > 1, i.e., the inequality Eq. (|4Tj) holds. If 
more than one phase satisfies this inequality, then we list 
all such phases regardless of the value of r\i . As discussed 
before, whenever two or more susceptibilities diverge, we 
cannot decide within our RG framework if the system 
chooses only one of these phases or if there is a coexis- 
tence. The resulting list of phases is shown in FigurelU In 
addition to the phases we found earlier in Section fill El 



namely the nematic (N), layer antiferromagnetic (AF), 
quantum spin Hall (QSH), layer-polarized (LP), Kekule 
current (KC), and s++ superconducting (s++ SC) states, 
a few other phases are predicted as possible outcomes of 
the RG flow if it ends on the target plane. These are: 
(a) Magnetoelectric phase (ME2): The order parame- 
ter for this phase transforms according to the E u charge 
representation. In this phase, currents forming a bow- 
tie pattern within a. plaqucttc appear. Like the nematic 
phase, this phase is gapless, but it reconstructs the low 
lying spectrum by lifting two of the four Dirac cones, (b) 
Kekule state (K): In this phase, a supercell made of three 
unit cells is formed, much like the Kekule current phase. 
The difference is that, in this phase, there are no currents. 
Instead, there is a modification of the hopping integrals 
such that the hoppings in one unit cell are unchanged, 
while, in the two other unit cells, the hoppings on alter- 
nating bonds are changed^. The phase is gapped, (c) 
Staggered spin current state (SSC): This phase is char- 
acterized by circulating spin currents in each plaquette 
flowing in opposite directions in each layer. This phase 
is not gapped, corresponds to a compensated semimetal, 
and the order parameter belongs to the A\ u spin repre- 
sentation, (d) s_| superconducting state (s.\ SC): Since 

a particle-particle susceptibility diverges in this case, a 
superconducting gap opens on both layers. The gaps 
are, however, not independent; they have opposite signs. 
The order parameter of this phase is a (charge 2) Az u 
spin singlet. 



Strictly speaking, when either x or y becomes infinite 
or they satisfy l + x-\-2y = 0, we are not allowed to divide 
by gE g as this coupling is not divergent. It shows up in 
Eq. (f73]) as a divergent A^E g )- Instead, these cases are 
explored by dividing by some other coupling. Wc follow 
the same procedure as described above in the case where 
we divided by gE g - Interestingly, since both A(E g ) and 
Bi,(E g ) diverge in the same way, the iji's are independent 
of the choice of the coupling that we divide by. 



In addition to the target plane, we also find the follow- 
ing four isolated stable fixed points. 
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with gE g (i — > oo ) > 0. In this case, only the ferro- 
magnetic (Ai g spin) susceptibility diverges. 
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and gA 2g (£ ~^ 00 ) < 0- The only divergent sus- 
ceptibility in this case is towards the anomalous 
quantum Hall stated {A2 9 charge). Here, charge 
currents circulate in each layer—, and in the same 
direction in both layers. 



<)., 



lim 



Vjy^A lu , 



(76) 



and gA lu (£ — > oo) < 0. This yields a loop current 
order—, or "orbital antiferromagnet" {A\ u charge). 
Like the above phase, there are charge currents 
circulating in each layer, but in opposite direc- 
tions. Note that the order parameter, T3I, can 
be thought of as a chemical potential shift with 
opposite signs in each valley. Therefore, at weak 
coupling, this phase corresponds to a compensated 
scmimetal with electron and hole pockets. 
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with gA lg (i — > 00) < 0. Although we would intu- 
itively expect this fixed point to favor a supercon- 
ducting state, we find no particle-particle suscepti- 
bilitcs diverging. Only the A\ g charge susceptibil- 
ity, or equivalcntly the electronic compressibility, 
diverges. Therefore, we conclude that the system 
enters a phase segregated state. 

We can now make a connection between the results 
obtained in previous sections and the analytic results ob- 
tained here. For the set of initial couplings and param- 
eters analyzed in Sees. MI DHlHFl the flow at t c always 
converges to the target plane, and never to any of the 
isolated points Rj . In the case of forward scattering only 
at £ = 0, and in the absence of trigonal warping, none of 
the couplings from the u and K representations are gen- 
erated. The flow always ends at the point, (x, y) = (0, 0), 
in the target plane, which corresponds to a pure nematic 
state. With trigonal warping included, u representation 



couplings are generated, even when we start with for- 
ward scattering only. However, we still do not generate 
any of the K representation couplings. This means that 
the end point of the flow at t c is restricted to the y = 
line in the target plane. Decreasing the initial coupling 
strength gA lg , while holding the bare V3 fixed, causes t c to 
decrease. At the same time, x increases. We always find 
that x < 1. As seen from Fig. HO these points correspond 
to a pure nematic order. However, as t c is lowered, the 
point in the target plane moves closer to x — 1, which is 
the intersection of the AF and QSH regions in the target 
plane. Upon reaching t c = exactly, we find that the ne- 
matic order parameter is absent while the AF and QSH 
susceptibilities become divergent. This is illustrated in 
Fig. EJb). Note that there is no point at which the only 
two diverging susceptibilities are AF and QSH, either in 
the target plane or as one of the four isolated points. This 
is because the asymptotic behavior of the $ functions is 
different at t = and our analysis, in which we assumed 
that t c > 0, does not apply there. 

When the RG flow begins with a finite backscattering, 
i.e., gE K 7^ 0, all 9 couplings are generated. We find that 
our numerical results correspond to points in the target 
plane where y ^ 0. With the physical constraints we 
impose on the initial couplings, |<?a 2 J < 9Ai„i lff£ K | < 
gA t /2 and gA ± positive, only the central region of the 
target plane is approached. With these constraints, we 
do not find any set of initial couplings for which K, SSC, 
or ME2 phases appear. 

In the previous sections, the flows never reach any of 
the isolated points Rj. However, one can expect that 
the flow will tend to one of these points if one starts 
with bare couplings sufficiently close to the associated 
ray and with a large initial value of I. To confirm this, 
we analyzed the flow equations with no trigonal warping 
and <7Ai = 0) < 0. For sufficiently large initial values 
of the interaction %g Aig {l = 0) < ^g c Mg « -0.13, 
the flow takes the couplings towards the R4 fixed ray. As 
stated before, this represents a compressibility instability. 
However, when g Alg (£ = 0) > g Al the couplings diverge 
toward the nematic fixed ratio, i.e., our flows end up on 
the target plane. 

The symmetry properties of the Hubbard limit have 
consequences for the asymptotic behavior of the RG flows 
studied in this section. When all three couplings, 5gi, 
Eqs. (|59"|) - (|6"Tj) are absent at the bare level, we have shown 
that they remain zero throughout the entire flow. As ar- 
gued above, the ratios of couplings at t c must lie either 
on the target plane or at one of the four isolated points 
Rj. In the target plane, the condition, 5gi = 0, is sat- 
isfied only when x = —2y. This defines a line of fixed 
points (strictly speaking a circle, since two points at in- 
finity are equivalent). As shown in Sec. IIIIF1 because 
Sgi = 0, the susceptibilities towards the layer-polarized 
and s++ superconducting states are identical. Therefore, 
this holds along the entire line x = —2y. For the repul- 
sive Hubbard model, a wide range of initial conditions, 
1CP 8 < t c < 1, maps onto the segment of this fixed line 
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FIG. 10: Plot of the critical temperature T c as a function 
of the layer energy difference V± , under such conditions that 
T c (0) = 0.01 A 2 /2m*. The critical temperature becomes zero 
when V± = Vj_,c = 4.93 x 10~ 3 A 2 /2m*, which corresponds to 
an applied electric field of ~ 16 mV/nm. 



that lies within the AF-only region. These results were 
also used in studying the attractive Hubbard model due 
to the U — > —U correspondence presented in Sec. IIII Fl 
The only difference is that both x and y change their 
sign under this mapping. The resulting fixed points are 
therefore situated in the region where the layer-polarized 
and s++ superconducting orders overlap. 

In addition to the fixed line that is part of the target 
plane, the condition that all three SgiS are zero is sat- 
isfied at the isolated fixed point Ry. However, we never 
find a flow toward that point for any set of bare couplings 
studied here. 



IV. EFFECT OF A PERPENDICULAR 
ELECTRIC FIELD ON THE PHASE 
BOUNDARIES 

We now consider the effect that applying a perpen- 
dicular electric field has on the phase boundaries of our 
system. This field creates an energy difference between 
the two layers of the sample, thus introducing a new term 
into the Hamiltonian, 



|k|<A<r=t4 



(78) 



We state the effects that this has on the Green's function 
and on the associated identities that we use in Appendix 
iBl and simply quote the main results here. The RG flow 
equations for the coupling constants become 

^=^9i9k^4?k*MQMQ,m, (79) 

j,k a—l 

where, in addition to the dimcnsionless parameters for 
the <f> functions that were defined before, we have one 
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The $ functions are given by Eqs. (|B3|) (|B15|) . 

In addition, the energy difference V± has a nontrivial 
behavior under rescaling; it obeys the flow equation, 



dv 



l+F(v 3 ,v,t)J2bi9i 



(81) 



where the coefficients 6j are given in Appendix C, and 
the function F is given by Eq. (|B17|) . 

We studied the behavior of the critical temperature 
as a function of v in a case where we know that the 
system enters the nematic phase when v = 0, namely 
when all coupling constants are zero except for gA lg > 0. 

A 2 

We assume that the energy scale ^— - = 200 meV, and 
that the critical temperature at zero field is T c = 2 meV; 
i.e. t c = 0.01. To determine t c for a given initial gA lg , 
we start with a high value of t and integrate our RG 
flow equations numerically up to a large value of £, say 
10. If we do not encounter any divergences in the flows, 
we lower t and integrate again. We continue until we 
find the highest temperature at which we encounter a 
divergence. Under the stated conditions, the behavior of 
t c as a function of v is as shown in Fig. 1101 

We find that, at v = 4.93 x 10 -3 , t c becomes zero. Con- 
verting this into an electric field using the stated energy 
scale and the formula relating the applied electric field 
to the size of the gap in the spectrum, A = 4^, where 
d is the distance between the two graphene layers and 
k w 3 is a factor accounting for imperfect screening^, 
we find that the electric field required to drive t c to zero 
is ~ 16 mV/nm. 

We have also determined which phase the system en- 
ters at all points on the curve in Fig. [TU]at which T c ^ 0. 
We did this by deriving the RG flow equations for the 
source terms and the formula for the free energy per unit 
area using the same procedures as before, but with prop- 
erly modified Green's functions, whose forms are given in 
Appendix B. We then numerically integrate the RG flow 
equations, and, from these solutions, determine the sus- 
ceptibilities to the phases corresponding to each source 
term just above the critical temperature; the phase with 
the highest susceptibility is considered to be the phase 
that is present. Using this procedure, we determined 
that, for all V± < Vj_ )C , the system enters the nematic 
phase. 



V. DISCUSSION 

The key result of this work is the identification of 
the conditions on the electron-electron interactions un- 
der which various electronic ordering tendencies, if any, 
dominate in half-filled bilaycr graphene. Our results for 
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the ordered states are summarized in Figs. [5] and [9j Aside 
from our use of one- loop RG, no further approximations 
are made. Therefore, our results can be stated rigor- 
ously at the level of mathematical theorems. While a 
large number of phases is, in principle, possible in the 
entire 9-dimensional space of couplings, as one can see 
from Fig. [HI our assertion is that the electronic nematic 
appears to be the unique dominant instability when for- 
ward scattering dominates. Similarly, the layer antiferro- 
magnet appears upon inclusion of sufficiently strong back 
and layer imbalance scattering. 

A similar approach in one spatial dimension^ results 
in divergences in the scattering amplitudes at finite tem- 
perature, naively suggesting a finite temperature phase 
transition, which we know cannot happen. Nevertheless, 
among many possibilities, the method does identify the 
correct channels for which long, but finite, correlation 
lengths develop. For example, the low-energy effective 
field theory for the course-grained half-filled Hubbard 
model does correctly determine that the dominant corre- 
lations appear in either the pairing (attractive U) or AF 
(repulsive U) channel^. Away from any special filling, a 
metallic state is also correctly predicted^. 

We view our RG results for the half-filled Hubbard 
model similarly. While the method correctly determines 
the dominant ordering tendency, there can be no finite 
temperature phase transition in 2D to either the AF state 
or the SC/LP state. A continuous spin SU(2) symmetry 
in the former case, or a continuous pseudospin symme- 
try in the latter case, would have to be broken at finite 
temperature, which we know cannot happen. Therefore, 
if the RG procedure had been performed exactly, none of 
these susceptibilities would have diverged as long as the 
temperature was finite. Interestingly, in this regard, the 
nematic state is different. This is because, when trigo- 
nal warping is included, as it is in our model, the broken 
rotational symmetry is discrete and thus it is possible 
to have a finite-temperature transition into this phase in 



2D. As argued previously-ii, this transition is continu- 
ous and belongs to the 3-state Potts model universality 
class^. Nevertheless, the non-mean-field exponents de- 
termined approximately from our fermionic model at one 
loop should not be expected to be accurate. It would be 
very interesting to see whether going to higher order ei- 
ther improves the accuracy of the exponents in the case 
of the nematic state or eliminates the finite-temperature 
phase transition altogether for the case of 0(3) order pa- 
rameters. The effects of (weak) disorder have not been 
addressed here either. These considerations may be im- 
portant in fully understanding the current experimental 
results^—. 

Even if the RG method used here is not without its 
limitations, it is unbiased and capable of systematically 
treating the leading instabilities in both particle-hole and 
pairing channels. In fact, for a large range of tempera- 
tures above the transition temperature, the couplings sat- 
urate to small finite values as all modes are eliminated, 
giving full justification to our method. In the special case 
in which our continuum field theory corresponds to the 
weak-coupling honeycomb bilaycr Hubbard model, we re- 
cover some of its non-trivial, exactly known, properties. 
This gives further support for the validity of our results. 
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Appendix A: Asymptotic behavior of the <E> functions 



The finite-temperature Green's function that is used to calculate the RG flows for infinitesimal symmetry breaking 
source terms is 



Gk(iw n ) 



iw„ls + ~ kl)lcril + v 3 k x T3(Jil + —k x k y T 3 a 2 l ~ vzk y \G 2 \ 



(Al) 



2 ti 



(1 + ST 3 ) 



IUJ,, 



«1 + (^zk 2 cos26> k + sw 3 fc cos 6^)0-! + (s^—k 2 sin26> k - u 3 fcsin 6» k )er 2 



v 2 k 2 + s^rv 3 k 3 cos 30 k 



■1. (A2) 



Throughout the Appendix, we will use the notation, TjO-jS^., for the (8 x 8) matrices that appear in our expressions; 
the Pauli matrices operate in valley, layer, and spin space, respectively. We find the following identity useful when 
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calculating the flow equations 

lA(l-di) 



f dkk 1 ^ d o k 

/ ~T~ r 1^ / ~o~ G k (zw„)(g)G ±k (±iw n ) 



8tt 



Tl 8 <8 1 8 [$i (>3, *) + $2 (i>3, *)] + 2 (^i 1 (8 lo-il + r 3 cr 2 l ® r 3 a 2 l) [$ 3 (^3, t) + $ 4 fa, *)] 

-T3I4 <g> T3I4 [$l (1/3, *) ~ $2 («3, *)] ± - (t 3 (Ti1 <g> T 3 (7il + lCT 2 l ® 1(T 2 1) [$4 (f3, *) - $3 (^3, *)] 



The $ functions arc defined as 



where 



*l(f3,*) = 

$2(^3, *) = 
$3^3, *) = 
$4(>3, *) = 

Ti(x, 1/3, t) 



ll/" 1 



2tt t 7-i Vl 3 ^ 
1 1 f 1 dx 1 



:Ti(x,i/ 3 ,t), 



T 2 (x,^3,i), 



11 — ^ 



2 /■! 



da; 1 



1 1 f 1 dx 



T 3 (x,vs,t), 



2tt t J_ t 



cosh 2 (%) 
T 2 (x, 1/3, = \Q\ tanh 

£1 qa v a ; 



T^tanhf^ 



A=± 



Ox 
2t 



T 3 (x,v 3 ,t) 
T 4 (x, 1/3, t) 



A=± 

-1 2i 

tanh 

cosh-M^ <9+ V 2 * 



(%) 



and 



Q 



1 + v\ ± 2xz/ 3 . 



In the limit of w 3 = we have 



$i(0, t) = $2(0, t) = tanh 
$3(0, i) = $4(0,0 = tanh 
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1 1 




2l + 
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2* cosh 2 


1 
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In the limit of T = we have 



$i(^3,0) = $4(^3,0) = 



4v 3 



$2(^3,0) 



TTl + Z/3 V(! + ^3) 2 
2^3 4^3 



2 1 (1- ^ 3 ) 2 

7T ^3 1 + ^3 V" V 1 + v l ' (! + ^3) 



n 



- # 



4^3 



(1 + 



.(A3) 

(A4) 
(A5) 
(A6) 
(A7) 

(A8) 

(A9) 
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(All) 



(A12) 

(A13) 
(A14) 

(A15) 
(A16) 
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Here, the complete elliptic integrals of the first, K(x), and the third, U(x,y), kind are defined as 



K(x) = I y (A18) 

o v 1 — x sin (p 



n(x,y) = / J (A19) 

J a (1 — ism 0)vl — 2/sm <p 

Physically, the logarithmic singularity associated with K(x) has its origin in the logarithmic divergence of the density 
of states at the van Hove point, where the lines of constant energy near each K point change from single to four closed 
contours. These log singularities appear only at t = 0, for t > they are smeared out. Because the divergences are 
intcgrable, they aren't the cause of divergence of g[s in (|T9l) . Instead, the coupling constants receive a "boost" at £ 
where v 3 {C) = 1- 



Appendix B: Green's functions in the presence of an applied perpendicular electric field 

In the presence of an applied electric field, the Green's function becomes 

G k (iuj n ) = (-iui n l 8 + d£lcril + v 3 k x T 3 oi\ + e^T 3 er 2 l - v 3 k v l<T 2 l + V±la 3 iy 1 
lr^, iu n l + (^fc 2 cos26> + sv 3 kcos8)a 1 + (s^fc 2 sin 26> - v 3 ksine)a 2 + V±1<J3 , > 

~ 2 hi w « + i^ fc4 + ^ + s^v 3 k* cos 39 + VI ' ( ' 



The generalization of Equation (|A3[) for this case is 



f kdk 1 ^ r*« d&k 

JA(i-M) 27r P„__Jo 2tt 



771 

— d£[(Tls (8 Is - T3I4 <E> t 3 1 4 )$i(^ 3 , v, t) + (=Fl 8 ® 1 8 + t 3 1 4 <8> nU)^ 2 {v 3 , v, t) 

+ i(lCTll ® lCTll =F 73(71 1 (gl T 3 fTll =F lC2l (8 l^l + T3tT 2 l ® T 3 2 l)<& 3 {v 3 ,V,t) 
+ \(l(T\l <g> lCTil ± T3CT1I (g) T 3 f7il ± lcr 2 l ® lcr 2 l + T3tT 2 l <8> T 3 er 2 l)$4(V 3 , t) 



where the $ functions are 



-(lcr 3 l ® I03I + T 3 (7 3 l (g> T3Cr 3 l)$ 5 (^3, U, t) 
+ (lcr 3 l ® lcr 3 l - r 3 cr 3 l (g) r 3 cr 3 l)*6(^3, M)L (B2) 
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ll/" 1 dx 
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The T functions arc 
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Y 3 (x,v 3 ,v,t) 

r 4 (x,v 3 ,v,t) 
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Q± 
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Q+ V 2t 
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u 2 ± 2xz/3, 



f| ± 2xv 3 . 



One other identity that we will find useful is 

kdk 1 
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2tt j3 
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2tv 



where F{v 3 , v, t) is 



F(v 3 ,v,t) 
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2tt 



■ tanh 



la3lF(iy 3 ,«,t) di, 



Q+ 
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(B12) 
(B13) 



(B14) 
(B15) 
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(B17) 



In arriving at this result, we will, in the intermediate steps, also find a term proportional to r 3 a 3 l. However, wc 
can use the periodicity of the integrand to show that this term will be zero at the end. We can also see that this 
must happen due to the symmetries of our system. Imagine that we tried calculating the expectation value of an 
observable, which would be represented by a matrix TiGjSk- This expectation value will only be non-zero if the matrix 
is proportional to one of the matrices appearing in the above identity, since said expectation value involves a trace 
of the product of the Green's function with the associated matrix. If a term proportional to T3CT3I were present, 
then this means that we can have a finite expectation value of the associated observable, which would, in this case, 
correspond to the gap opened by an anomalous quantum Hall order parameter. This order parameter breaks time 
reversal symmetry. However, we should not be able to develop a finite expectation value of this observable because 
our Hamiltonian is time reversal invariant. ^ 



Appendix C: Determination of RG flow equations 

We now show how to derive the flow equations for the 
four-fermion coupling constants g and the source term 
constants A. We start by performing a cumulant ex- 
pansion of the partition function to second order in the 
"perturbation" Sj nt + AS. 

Z » cxp [-(S int + AS) Q + i((5 mt + AS) 2 ) 0iC ] , (CI) 

where, in the subscripts on the averages (•), "0" means 
to average with respect to the bare action So, and "C" 



means that the average is "connected" ; that is, it can be 
represented with connected Feynman diagrams. We now 
integrate out modes in thin shells; by doing so, we gen- 
erate terms that renormalize different constants in our 
theory. Wc will first discuss the terms that renormal- 
ize the four-fermion coupling constants, since the general 
procedure is the same. There are five different types of 
terms that appear; these are represented by the diagrams 
shown in Figure 1111 
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FIG. 11: Diagrams representing contributions to the renor- 
malization of the four-fermion coupling constants (?;. The 
dashed lines represent 8x8 matrices, the black lines repre- 
sent slow modes, and the red lines represent fast modes. 




FIG. 12: (a) Diagrams representing contributions to the renormalization of the particle-hole source terms. All lines are as in 
Figure[TT] In addition, the wavy lines represent the source terms, (b) Diagram representing contributions to the renormalization 
of the particle-particle source terms. 



The first diagram gives the following correction: 



/ \l Tr[51G k (io;)f/lG k (ia;)] 1^^(1)^(2)^(3)^(4), 
s.u Jl - 2 ' 3 - 4 > <r,*> 



(C2) 



where the numbers 1 — 4 are shorthand for the momentum and frequency dependences of the Grassman fields, and 
/i 2 3 4 re P resen ts the integrals and sums over these variables along with the proper momentum- and frequency- 
conserving 8 functions, and similarly for J k bj . We may evaluate the integral and sum over k and u, respectively, 
using Eq. (|A3[) in the absence of an external electric field or Eq. (|B2[) when said field is present. In both cases, 
this term is only nonzero if S = U, so that we generate a correction to gs that is proportional to g 2 s . The nonzero 
contributions to the coefficients, 



41 = A §(1) + 4S( 2 + 3 ) + A S( 4 ) + A S(5), (C3) 



in Eq. (Q35J) are 



4/ 2) (i) = -Hsii^rfVsU) 2 ]}^, (C4) 

= i{Tr[(rf harl) 2 ] T Vil) 2 ] T Tt[(lf >la 2 l) 2 ] + Tr[(lf V 3 <7 2 1) 2 ]}^, (C5) 

4 /6) (l) = |{Tr[(rfh<73l) 2 ]±^[(rfV 3f r3l) 2 ]}J. (C6) 

In these expressions, the top signs correspond to the first number in the superscript on the left-hand side, while the 
bottom corresponds to the second. The coefficients only enter into our analysis when a finite electric field is 
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present. The 8x8 matrices are defined as follows: 



rr = is (C7) 

lf> - r 3 a 3 l (C8) 

r« = lail, T 3 2) =r 3 a 2 l (C9) 

ri 1} = r 3 i 4 (cio) 

rl 1 ' = ltr 3 l (Cll) 

lf } = r 3C 7il, r^ 2) = -la 2 l (C12) 

4 1} = nail, 4 2) =T 2 ail (C13) 

if } = no-al, 4 2) =r 2( r 2 l (C14) 

if > = nu, if 5 = -r 2 a 3 i, if 5 = -r 2 i 4 , if 5 = - nt T 3 i. (Ci5) 

The superscripts (m) refer to the multiplicity of a given representation. Here, and throughout this appendix, A^^{ri) 

represents the contribution to from the nth diagram in Fig. QT] 
The second and third diagrams together give the following correction: 



5S 2+3 = -Y,9s9u f ^(1)5^(2)^(3) f UlG k (iu)SlG k (iuj)Ul 

su ^1,2,3,4 \J\s.^.uj 



V(4). (C16) 



Note that the first two ip fields carry an explicit spin index. The second two do not; for notational simplicity, these two 
are extended to be eight-component spinors in valley, layer, and spin space. In both cases that we consider, the second 
matrix U\G k {iui)SlG k (iu3)Ul appearing in this expression is proportional to 51. Therefore, this term also represents 
a correction to gs, but now it generates terms involving the products, gs9u- We may extract the contributions to 
the A^l coefficients by noting that, since the second matrix is proportional to 51. Using Tr(r i - m ' ) r^™' > ) = 8SijS mn , we 

find that the nonzero contributions to the A^' k coefficients are 

m i * 

4/ 2) (2 + 3) = \ £{Tr[(rWrf >) 2 ] ± Tr(rf>rf Wf V 3 l 4 rf>)}|-, (C17) 

rn—1 

rrij 

Ag/ 4) (2 + 3) = -i ^[TrCrWrfhaxirWi^irf^T^^rf ViiifViiif*') 

m— 1 

T Tr^^r^l^ir^l^irf 1 ') + Tr(lf V 3 a 2 irfV 3 a 2 irf >)]^, (C18) 
A% 6) (2 + 3) = [Tr(lf ^fh^irf Wlf°) ± Tr(r«rf ) T3 a 3 irf V 3( 7 3 ir( m) )]^-. (C19) 

m— 1 

Here, the sum on m is taken over the multiplicity of the jth representation, and the "2 + 3" in the "arguments" means 
that the given contribution is the total contribution from the second and third diagrams. 
Finally, the fourth and fifth diagrams give the following: 

SS 4 = -\Y,9s9u I f ^(l)SG k (iu)U^j(2)^(3)UG k (iuj)S^(4) (C20) 

§ U * 1)2,3,4 «/k> ,w 

6S 5 = [ f ^{l)SG k {iu)U^{2)^{i)SG^{-iuj)U^) (C21) 

SU Jl,2,3,4 Jk>,w 

Both matrices occurring in each expression are proportional to each other, but will in general not be proportional to 
cither 5 or U . These terms therefore represent corrections to a coupling gy that are proportional to gsgu- Using the 
same observation as before, we can find the contributions to the A^l coefficients. Denoting V — Tu, these are 

mi raj 

, m 



< /2) (4) = E E[ Tr ( r i 1)r ^ )r 5" ) ) Tl -( r i 1)r 5" )r ' m) ) ± Tr(rWr|™V 3 i 4 r("))Tr(r«r("V 3 i 4 r( m ))]^, (C22) 

rn—1 n— 1 

m l rnj 

4l /4) ( 4 ) = -ale E E^rj^r^Wr^)^^ 



rn—1 n—1 
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t Tr(r^ m ha 2 irf l) )Tr(r^r'fha 2 irf n) ) + TV(r«r^V 3CT2 irf ^(r^rf Vaasir.H)]^-, (C23) 

m— 1 n— 1 

(C24) 



and 
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m—l n—1 

mi rrij 

^i 3/4) ( 5 ) _ V^rTV/T>( 1 M m )-i~ir>( T, hl2 _u mrwrWrt 1 ™),^^ 1 r(")M2 



256 / 
m—l n—1 

[Tl . (r W r (™) 1( 

mi nij 



p(i)p( m ; 




i 


) T 3 (Tl 


f)] 2 + 


[Tr(r (i) r ("0 T3(J2ir (n: 


)] 2 } 


m* 

47r' 




Wr< n) )] 2 ±[Tr(r< r 


lTi(m 

i 


V3C73 



m 
4-7T 



(C26) 
(C27) 



We now turn our attention to the symmetry-breaking source terms. In this case, we have different procedures for 
the case without an applied electric field and the case with one. We will consider the former case first. The corrections 
to the particle-hole and particle-particle source terms are represented by the diagrams in Fig. 1121 

The particle-hole source term corrections give us 

SS vh = EE A ^-9s / / Tr[G k ( lW )0«G k (^)51]^(w')^lV'k'(w') 

- EE A "^ / / i>l("')SlG k (iuj)O^G k (iu))SlTPv(u)'). (C28) 

In the first term, the trace will only be nonzero if £1 = 0^ % \ and, in the second term, the matrix appearing in the 
expression is proportional to O' 1 ' . Therefore we see that different source terms are not mixed to this order. Note that 
the first term is only nonzero if represents a charge order, and vanishes for spin orders. The contributions to the 
coefficients 

BW=flW(l) + flW(2) (C29) 

in Equation (|28|l are 

rrij 

B^ /2) {1) = -^MOWrf^Tr^l^rsUrf)]^, (C30) 



47T 

Bf/^iX) = iEtT^^i^^^i^^Tr^aaoWraaiir^) 

n=l 



T Tr(la 2 10«l ( r 2 ir( n) ) + Tr(r 3( 7 2 10«r3 ( 7 2 ir5" ) )]^-, (C31) 
i?|; /2) (2) = ig{Tr[(0( i )r(" ) ) 2 ]±Tr(OWrf ) r3l40 ( ' i) r 3 l4rf ) )}^, (C32) 



n=l 
m 

4 3/4) (2) = E[^(C (i) rf haxlO^laxlY^) T Tr(0»rj" ) T 3 a 1 10«T3 CT iirf >) 

n=l 

T Tt(0«rf Wo (i Wlf°) + Tr(0«rf Va^lOW^^ir^)]^. (C33) 



Here, the "arguments" have the same meaning as before, but with respect to Fig. [T2l 
The correction to the particle-particle source term is 



«V = -|EE A ^ / / ^("')^GkMO (i) [G- k H W f(51) r C k -(^')+cc. (C34) 
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For similar reasons as above, the product of five matrices appearing in this expression is proportional to OW ; and 
therefore different source terms are not mixed to this order. Also note that the 8x8 matrix OW must be completely 
antisymmetric. The values of the coefficients B^' in Equation (|29|) are therefore 

B% /2) = -^^{Tr^WrfdW^f^J^TrtOWrfVa^OWra^Cr^f]}^-, (C35) 

Bg /4) = -i ^{Tr[0( i; )r(" ) la 1 10( l )l ( T 1 l(r;" ) ) T ] ± Tr[0«r;" ) r 3 a 1 10( i; )r 3 ail(r; rl) ) T ] 

n=l 

T TrlOWr^l^lOWl^l^) 7 ] _Tr[OWr(" ) r 3 a 2 ld(' i )r 3( T 2 l(r(" ) ) T ]}^. (C36) 



Now we consider corrections to the finite applied elec- 
tric field; see Eq. (|5Tj) . In this case, we find that the 
lowest-order corrections come from the first-order term 
in the cumulant expansion; these first-order corrections 
would be zero in the absence of the electric field. They 
are represented by "tadpole" and "sunrise" diagrams, as 
shown in Figure[T3J The contribution from the "tadpole" 
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FIG. 13: Diagrams representing contributions to the renor- 
malization of the applied electric field term. All lines are as 
in Figure [12] 

diagrams is 

6S t = -J2as [ f Tr[SlG k (ia;)]4,(a/)SlV>k<(w')- 

(C37) 

The integral over k and sum over u> can be evaluated 
using Equation (|B16I) . The trace occurring in this ex- 
pression is only nonzero if S = las. Therefore, we only 
generate a correction to the applied electric field. Since 
ler 3 l = we see that the only non-zero contribution 
from this term to the coefficients bi in Eq. (|8Tj) is to 65, 
and this contribution is 65 (tadpole) = 8 x 
The "sunrise" diagrams give us 

5S S = Y,9S I I 4V)<SlGkMSlVv(c<A 

g ./k^.cj' J]s. > ,uj 

(C38) 

The matrix occurring in this expression is proportional to 
ler 3 l, and thus we, once again, only generate corrections 
to the applied electric field. This will contribute to all of 
the hi. These contributions are given by 

^(sunrise) = \ £ ^{\a^ n) I^IY^)^. (C39) 

m 



The total value of bi is simply the sum of the above two 
contributions, i.e., bi = bi (tadpole) + bi (sunrise). 

Appendix D: Coefficients in the free energy 

The coefficients appearing in the free energy, Eq. 
(|37)) . are 

af /2i = 8±Tr[(0«r 3 l 4 ) 2 ], (Dl) 

<4%,i = -H T r[(0 W l ( T 1 l) 2 ]TTr[(0«r 3 a 1 l) 2 ] 

T Tr[(0(*W) 2 ]+Tr[(0«T 3 a 2 l) 2 ]}. (D2) 
The c?^ { coefficients are 

</2,i = 8TTr[(0«r 3 l 4 ) 2 ], (D3) 

T Tr[(0«la 2 l) 2 ]-Tr[(0«r 3( 7 2 l) 2 ]}. (D4) 

Appendix E: Analytic determination of the phase 
boundaries in the fixed ratio plane 

We will now describe how we determined the phase 
boundaries in the target plane. These boundaries are 
defined by the sign of the susceptibility exponent ji , as 
given by Equation (|5tJ)) , for a given phase; whenever it 
is positive, we say that the associated phase is present. 
The value of A(E g ) is given by Equation (f73|) . We may 
obtain B^ie ) from Eqs. (|35p and (|3l))) and from the 

coupling constant ratios p\ 9 given in Equations (|70|) - 
(|72l . Because of this, all of the ji will have the form, 

= Qi(x,y) , . 

'* 3 + 2x + 3x 2 + 4y + Axy + 8y 2 ' 1 ' 

where Qi(x, y) is an inhomogencous quadratic function of 
x and y. The denominator of this expression is positive 
definite, so that the sign of the exponent is determined 
entirely by Qi(x,y). Our condition that 7$ be positive 
thus requires that Qi(x,y) > 0. We therefore see that 
the phase boundaries, given by Qi(x, y) = 0, are all conic 
sections. 
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